

%% "true" DGP point estimates
load('../true_DGP/tDGP.mat','DGP')
true_beta=DGP(1:(2*N));
clear DGP


%% reset rng seed
rng('default')


%% Great Depression, 1929-1933
data=readtable('../../Data/data_main_sample.csv');
deep_shock=mean(data.y(data.year>=1931&data.year<=1931&isfinite(data.y)));
shallow_shock=mean(data.y(data.year>=1929&data.year<=1933&isfinite(data.y)));
clear data


%%
NSP=1e+3; % number of simulated sample paths
epsilon=mvnrnd(zeros(N,1),Sigma,10*NSP)';
epsilon=reshape(epsilon,N,10,NSP); % growth shocks
pD_gd=zeros(NSP,10);


%% short-deep depression
T_gd=1;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    y_gd(:,1)=deep_shock+ep(:,1); %#ok<*PFBNS>
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            y_gd(:,t)=deep_shock+ep(:,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_sd1=mean(pD_gd);


%% long-shallow
T_gd=5;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    y_gd(:,1)=shallow_shock+ep(:,1); %#ok<*PFBNS>
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            y_gd(:,t)=shallow_shock+ep(:,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_ls1=mean(pD_gd);


%% short-deep, concentrated among autocracies
T_gd=1;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    shock=deep_shock*sum(M(:,50)==1)/(2*sum(D(M(:,50)==1,50)==0)+sum(D(M(:,50)==1,50)==1));
    y_gd(D(:,50)==0,1)=2*shock+ep(D(:,50)==0,1);
    y_gd(D(:,50)==1,1)=shock+ep(D(:,50)==1,1);
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            shock=deep_shock*sum(M(:,50+t-1)==1)/...
                (2*sum(D_gd(M(:,50+t-1)==1,t-1)==0)+sum(D_gd(M(:,50+t-1)==1,t-1)==1));
            y_gd(D_gd(:,t-1)==0,t)=2*shock+ep(D_gd(:,t-1)==0,t);
            y_gd(D_gd(:,t-1)==1,t)=shock+ep(D_gd(:,t-1)==1,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_sd05=mean(pD_gd);


%% long-shallow, concentrated among autocracies
T_gd=5;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    shock=shallow_shock*sum(M(:,50)==1)/(2*sum(D(M(:,50)==1,50)==0)+sum(D(M(:,50)==1,50)==1));
    y_gd(D(:,50)==0,1)=2*shock+ep(D(:,50)==0,1);
    y_gd(D(:,50)==1,1)=shock+ep(D(:,50)==1,1);
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            shock=shallow_shock*sum(M(:,50+t-1)==1)/...
                (2*sum(D_gd(M(:,50+t-1)==1,t-1)==0)+sum(D_gd(M(:,50+t-1)==1,t-1)==1));
            y_gd(D_gd(:,t-1)==0,t)=2*shock+ep(D_gd(:,t-1)==0,t);
            y_gd(D_gd(:,t-1)==1,t)=shock+ep(D_gd(:,t-1)==1,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_ls05=mean(pD_gd);


%% short-deep, concentrated among democracies
T_gd=1;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    shock=deep_shock*sum(M(:,50)==1)/(2*sum(D(M(:,50)==1,50)==1)+sum(D(M(:,50)==1,50)==0));
    y_gd(D(:,50)==0,1)=shock+ep(D(:,50)==0,1);
    y_gd(D(:,50)==1,1)=2*shock+ep(D(:,50)==1,1);
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            shock=deep_shock*sum(M(:,50+t-1)==1)/...
                (2*sum(D_gd(M(:,50+t-1)==1,t-1)==1)+sum(D_gd(M(:,50+t-1)==1,t-1)==0));
            y_gd(D_gd(:,t-1)==0,t)=shock+ep(D_gd(:,t-1)==0,t);
            y_gd(D_gd(:,t-1)==1,t)=2*shock+ep(D_gd(:,t-1)==1,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_sd2=mean(pD_gd);


%% long-shallow, concentrated among democracies
T_gd=5;
parfor nsp=1:NSP
    ep=epsilon(:,:,nsp);
    y_gd=zeros(N,10); % 2000-2009
    shock=shallow_shock*sum(M(:,50)==1)/(2*sum(D(M(:,50)==1,50)==1)+sum(D(M(:,50)==1,50)==0));
    y_gd(D(:,50)==0,1)=shock+ep(D(:,50)==0,1);
    y_gd(D(:,50)==1,1)=2*shock+ep(D(:,50)==1,1);
    Pd_gd=zeros(2*N,10);
    B_gd=zeros(2*N,10);
    P_gd=P(:,:,50);
    % 2001
    nm=find(M(:,50)==1);
    DD=[diag(1-D(nm,50)),diag(D(nm,50))];
    P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
    Pd_gd(:,1)=sqrt(diag(P_gd\eye(2*N)));
    B_gd(:,1)=B(:,50);
    B_gd([nm;N+nm],1)=B_gd([nm;N+nm],1)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,1)-DD*B_gd([nm;N+nm],1))));
    U0=exp(alpha+theta(1)*(Pd_gd(1:N,1)*sgi_nodes(:,1)'+B_gd(1:N,1)+ep_sd*sgi_nodes(:,2)'));
    U0=(U0./(1+U0))*sgi_weights;
    U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),1)*sgi_nodes(:,1)'+B_gd(N+(1:N),1)+ep_sd*sgi_nodes(:,2)')-f);
    U1=(U1./(1+U1))*sgi_weights;
    D_gd=zeros(N,10);
    D_gd(:,1)=M(:,51).*(U1>=U0);
    % 2002-2010
    for t=2:10
        if t<=T_gd
            shock=shallow_shock*sum(M(:,50+t-1)==1)/...
                (2*sum(D_gd(M(:,50+t-1)==1,t-1)==1)+sum(D_gd(M(:,50+t-1)==1,t-1)==0));
            y_gd(D_gd(:,t-1)==0,t)=shock+ep(D_gd(:,t-1)==0,t);
            y_gd(D_gd(:,t-1)==1,t)=2*shock+ep(D_gd(:,t-1)==1,t);
        else
            y_gd(:,t)=true_beta(1:N).*(1-D_gd(:,t-1))+true_beta(N+(1:N)).*D_gd(:,t-1)+ep(:,t);
            nm=find(D(:,50+t-1)==D_gd(:,t-1));
            y_gd(nm,t)=y(nm,50+t-1);
        end
        nm=find(M(:,50+t-1)==1);
        DD=[diag(1-D_gd(nm,t-1)),diag(D_gd(nm,t-1))];
        P_gd([nm;N+nm],[nm;N+nm])=P_gd([nm;N+nm],[nm;N+nm])+DD'*(Sigma(nm,nm)\DD);
        Pd_gd(:,t)=sqrt(diag(P_gd\eye(2*N)));
        B_gd(:,t)=B_gd(:,t-1);
        B_gd([nm;N+nm],t)=B_gd([nm;N+nm],t)+P_gd([nm;N+nm],[nm;N+nm])\(DD'*(Sigma(nm,nm)\(y_gd(nm,t)-DD*B_gd([nm;N+nm],t))));
        U0=exp(alpha+theta(1)*(Pd_gd(1:N,t)*sgi_nodes(:,1)'+B_gd(1:N,t)+ep_sd*sgi_nodes(:,2)'));
        U0=(U0./(1+U0))*sgi_weights;
        U1=exp(alpha+theta(2)*(Pd_gd(N+(1:N),t)*sgi_nodes(:,1)'+B_gd(N+(1:N),t)+ep_sd*sgi_nodes(:,2)')-f);
        U1=(U1./(1+U1))*sgi_weights;
        D_gd(:,t)=M(:,50+t).*(U1>=U0);
    end
    pD_gd(nsp,:)=sum(D_gd)./sum(M(:,51:end));
end
pD_gd_ls2=mean(pD_gd);


%% for Figure 6
results=table(years(51:60),pD(51:60),pD_hat(51:60),pD_gd_sd1',pD_gd_ls1',...
    pD_gd_sd05',pD_gd_ls05',pD_gd_sd2',pD_gd_ls2',...
    'VariableNames',{'year','data','model','sd','ls',...
    'sd_autoc','ls_autoc','sd_democ','ls_democ'});
writetable(results,'../../Figures/Figure_6/great_depression.csv')



